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Nonlinear System Guidance in the Presence of Transmission 

Zero Dynamics 


G. Meyer, L. R. Hunt, and R. Su 
Ames Research Center 

Summary 


An iterative procedure is proposed for computing the commanded state trajectories and controls 
that guide a possibly multiaxis, time-varying, nonlinear system with transmission zero dynamics through 
a given arbitrary sequence of control points. The procedure is initialized by the system inverse with the 
transmission zero effects nulled out. Then the “steady state” solution of the perturbation model with the 
transmission zero dynamics intact is computed and used to correct the initial zero-free solution. Both 
time domain and frequency domain methods are presented for computing the steady state solutions of 

the possibly nonminimum phase transmission zero dynamics. The procedure is illustrated by means of 
linear and nonlinear examples. 


1 Introduction 

This report presents a procedure for guiding a possibly nonlinear system through a schedule of 
control points. The paradigm is a fully automatic aircraft subject to air traffic control (ATC). The 
ATC provides a sequence of waypoints through which the aircraft trajectory must pass. The waypoints 
typically specify time, position, and velocity. The time separation between the waypoints is normally 
greater than one minute. The flight vehicle management system (FVMS) on board the aircraft is 
normally aware of at least one waypoint in advance. The FVMS planner provides a sequence of control 
points between the waypoints so that the aircraft will satisfy additional, vehicle specific constraints. The 
guidance command generator transforms the schedule of control points into reference state trajectory 
segments that are flyable and that pass through the control points. These reference trajectories are 
solutions to the system state equation. The rest of the system is a model follower in which the regulator 
maintains the aircraft state close to the reference state. Thus, the structure of the complete system is 
hierarchical: at the top is the ATC; at the bottom are actuators. There is a progressive filling in of the 
detail as one moves down the hierarchy. In the present report we are concerned with algorithms for 
transforming the control point schedule into reference state trajectories. 

If the system to be controlled is approximately in pure feedback form, that is, if the zero dynamics 
are negligible, then system inversion provides an effective procedure for the generation of guidance 
commands (refs. 1 and 2). In their prize paper (ref. 3), Isidori and Byrnes solve the guidance problem 
m the presence of zero dynamics. However, a solution to a nonlinear partial differential equation must 
be obtained, and that is not always practical. Furthermore, the commanded output is restricted to outputs 
of a relaxing autonomous exosystem, whereas we are interested in the case where the output may have 
discontinuous higher derivatives at the control points. Paden, Chen, and Devasia (refs 4-6) made a 


major advance by finding an iterative solution that avoids the partial differential equation and admits 
discontinuities at control points. 


The solution described in the present paper may be viewed as a modification of the iteratio 
reference 6. Our procedure consists of two nested iterations. The computation is initialized by he pure 
feedback solution in which the zero dynamics are nulled out. Then the error resulting from the zero 
dynamics is removed (outer iteration) by solving a sequence of linear problems. Each linear problem 
is solved by means of the linear version of the iteration in reference S. The introduction o e ou 
iteration allows one to control the effects of nonlinearities on the convergence of the inner iteratiom 
The effectiveness of the approach is illustrated by means of applications to linear, time-varying, and 


nonlinear systems. 
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2 System Guidance 


The purpose of the guidance subsystem is to generate reference trajectories that will comply with 
the specifications of a higher-level planner and with the capabilities of the controlled process. Thus, 
for example, consider a simplified model of the interaction between the ATC and the aircraft under 
its influence. For simplicity, let both the ATC and all the aircraft be fully automatic. A model of the 
combined process is shown in figure 1. Based on the available airspace and goals and capabilities of the 
various aircraft, the ATC provides to each aircraft a sequence of waypoints. The FVMS of each aircraft 
develops a plan for passing through the assigned waypoints. Then the planner provides a control point 
table (CPT) of control points that specify the desired conditions at particular times, as well as the aircraft 
configuration in which the segments are to be flown. The command synthesizer connects the control 
points by functions of time from a standard set, such as polynomials. The result is the commanded 
motion of the output y c and the required time variation of the configuration variables p, such as flaps and 
landing gear. The guidance generator transforms the desired motion and configuration change ( UoP ) 
into the complete guidance state trajectory and control ( x 9, u 9). The plant regulator transforms the 
tracking error x p — x 5 into a corrective action u e , which is combined with u s to produce the plant 
control u p . The primary virtue of this model-follower structure is that the regulator is responsible only 
for controlling the uncertainty arising from disturbance and modeling errors, and not for shaping the 
plant response. It may be noted that the overall structure is hierarchical with regard to the horizon 
width and plan refinement and stability. ATC is at the top with the widest horizon, slowest-changing 
plan, and coarsest specification. The FVMS planner has a narrower horizon and a faster-changing, more 
refined plan. The guidance generator fills in a lot of detail, and the regulator is at the bottom’: it is 
purely reactive, with the narrowest horizon, and it produces the most refined commands. The ATC plan 
is stable for minutes, the guidance for seconds, and the regulator for milliseconds, changing at every 
sample (e.g., 20 msec). In the present report, we are concerned with design of the guidance generator, 
that is, with the computation of ( x 9 ,u g ). 


A high-level description of the guidance problem is straightforward. We are given the system state 
equation and the output map. 


x = f(x,u,p ) 
y = h(x,p ) 


( 2 - 1 ) 


where the state x £ R n , the control u £ R m , the parameter p £ R k , and the output y £ R m . We 



Figure 1. A model of the interaction between ATC and aircraft. 
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are also given the evolution of the parameter {p(t),t G T) and the desired evolution of the output 
{y c {t),t G T}. The problem is to find a control {u 9 {t),t € T} and the corresponding solution 

{x 9 (t), t G 7"} of the state equation so that 

h[x 9 (t),p(t)] = y c {t) (2- 2 ) 

Furthermore, the solution {x 9 ,u 9 ) should not reach values that are excessive for the task at hand, 
and the guidance state x 9 must be everywhere continuous. It may be noted that since the initial condition 
x 9 0 is not given, the guidance problem is not explicitly an initial value problem. 

The following very simple example illustrates the salient features of the guidance problem and the 
solution approach that will be pursued in the present report. 

Example 2-1 

Consider a scalar second-order system, 


X\ = X2 

±2 = u (2-3) 

y = -xi + X2 


that is supposed to track 


y c (t ) = sin ut 


(2-4) 


Obtain a differential equation describing the zero dynamics. Note that 
there is a nonminimum phase zero at s = 1. 

z — z = sinaif 
x i = z 
x 2 = z 
u — z 

Consider the initial value problem with the initial condition to be 
selected later, 

z — z — sin cut 
z(0) = ZQ 


The solution consists of the sum of the homogeneous and particular 


parts: 


Z = z h + z p 


(2-5) 


( 2 - 6 ) 


(2-7) 


where / t 

Zhtt) = ce ( 2 - 8 ) 

z p (t) = asin(t + 0) 
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and where (a, 9) are the 
zero dynamics. 

gain and phase of the transfer function of the 

aei® = ( ju > — l)' 1 

(2-9) 

Consequently, 

z(t) = ce t +a sin(i -f 9) 

(2-10) 


We choose z( 0) — z p ( 0) in order to obtain c = 0 and thus null the 
homogeneous part. The guidance trajectory is then given by the par- 
ticular solution, 

x\ p (t) = asin(cuf + 9) 

X 2 p (t) = a<jJcos(ut + 9) 
u p (t) = — acu 2 sin(c ot + 9) 


( 2 - 11 ) 


The following points should be emphasized: 

• In general, the reference y c as well as the guidance solution may persist indefinitely. 

.The reference y c can be expected to be differentiable at the control points only a small number 
or times, but the guidance state should be continuous everywhere. 

.In the constant linear case, the guidance solution is the particular solution that can be obtained 
by means of the (possibly unstable) transfer function of the zero dynamics. 

• We wish to avoid the computation of the guidance solution by numerical integration of the 
(possibly unstable) zero dynamics. 

• Whereas regulation algorithms must be causal, guidance, being a planning activity, can have 
noncausal computations. 


In the present report we generalize this approach to time-varying and nonlinear systems. The 

specific case of aircraft guidance is considered, but it should be clear that the methods are applicable 
to a wider class of systems. 
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3 Aircraft Model 


A rigid body model of an aircraft can be described by the following set of equations: 

t T — v r 

Vr = Clf[ ( u m , u p , u f , vf,u brb ) + g6 3 

C b r ^ (u> brb ')C br 


( 3 - 1 ) 


^brb ~ I™ i. u m ,U p ,Uf,vf, U> brb ) 

where r r ,v T G R 3 are the runway coordinates of the aircraft center of mass position and velocity, 
respectively; C br is the direction cosine matrix locating body-fixed axes with respect to the runway; 
Lo br b are the b°dy coordinates of the aircraft angular velocity relative to the runway; u m g R 3 is the 
moment control; u p G R 3 controls the magnitude and direction of engine thrust; and ut G R 3 represents 
the aircraft configuration variables such as flaps. Normally not all coordinates of (u p , ut) are available 
for control. The partition into active controls and configuration variables, respectively, is defined by the 
control mode. Finally, g is the acceleration of gravity, vf = C br (v r - w r ) are the body coordinates of 
the air velocity, and w r are runway coordinates of wind. The force and moment characteristics of the 
aircraft are described by the force and moment functions (//,/£"). They are typically multiaxis, highly 
coupled, and nonlinear. It should be noted that the state space is not flat, so coordinate patching may be 
required for some aircraft maneuvers. A convenient coordinate system is obtained by the introduction 
of an intermediate axis system, such as the nominal stability axis, so that 


— C b tCtr ( 3 - 2 ) 

where C bt is parameterized by Euler angles a bt in a given sequence, and C tr is an explicit function of 
time specifying the nominal attitude of the aircraft. The state equation in this coordinate patch becomes 

T 'jr — Vf 


v r — Cj T ff{u m , Up, uj, a bt , vf, a bt , t) + g6 3 


( a bt ) = 


( 3 - 3 ) 


(®bt) — ®bt — f a { u mi Up, uj, a b t,vf, a b i, t ) 

The state space is flat. The force and moment functions depend explicitly on time by virtue 

of the explicit time dependence of C tr . Three derived functions are useful for control system design 

purposes. The first is the moment trim map, which will be defined as a partial inverse of f a and 
denoted by g a : 

u m = 9 a (a bv up, u f , a bt , vf, a bt , t ) 

( 3 - 4 ) 

®bt ~ f ( w m> u Pi u fi ^bti v t > ®bti 0 
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so that g a computes the value u c m of the moment control that is required for the generation of a given 
Euler angle acceleration a c bt , whereas all the other relevant variables are held fixed at some given values. 
The second is the force function subject to moment trim. It will be denoted by F and defined by the 
composite 

F f (a c bt ,Up,u f ,a M ,vf,au,t) = f , [g a (a bt ,u p ,Uf,a bt ,vf,a bt ,t),Up,Uf,a bt ,Vt,a u ,t] (3-5) 
The third derived function, denoted as F°. is the force trim subject to moment trim. Here it will be 
defined implicitly as the solution to 

Ff{a c bt ,u p ,Uf,a bt ,Vt,a bt ,t) - ft = 0 


The control mode defines how {u p ,u f ,a bt ) is to be partitioned into dependent and independent 
variables For example, in the mode in which force is trimmed by means of roll, pitch, and thrust the 
function F a outputs the values of these variables for the given force ff and all the remaining variables. 


It may be noted that the selection of the coordinatization of the underlying state space and the 
selection of the control mode has a major effect on the form of the state equation. In addition, the 

output map h, 

y = h(rr,V r ,C br ,U brb> Up,Uf,W r ,C r t) v ' > 

is also defined by a mode-the output mode. Thus, for example, the output may be position f>. or 
velocity Ur or its spherical or polar coordinates; or it may be the relative air velocity u t , or perhaps 
aircraft attitude, C,„, or something else. Many output modes are of practical interest. The three ypes 
of modes namely the mode that defines the coordinate patch, the mode that defines the flow of the 
control action through the system, and the mode that defines the output function, will be considered 
components of a three-dimensional mode, which will be called the operation mode. Clearly, there are 
many values of the operation mode, and for each such mode there is a local state space model of the 
controlled process. The (possibly) multiaxis (multi-input multi-output) models have the following form. 

±1 = /l(zi,X2>P) 

X 3 = X 4 ^ ‘ 

XA = /4(xi,X2,X3,X4,U,p) 


y = xi 

where x h u € R m , and p € R k is a parameter, which is in general a given function of time. The state 
space is flat. The variables have been sorted out as follows. 


(1) is invertible with respect to X 2 , 

x 2 = 9i(xi,±i,p) _ (3-9) 

fi[xi,gi(x\,±i,p),p] = xi 


(2) U is invertible with respect to the moment control variable u. The moment trim map g in 
equation (3-4) will be denoted by 

u = <74(xi,X2,X3,X4,U2iP) (3-10) 

f A [x 1 , X 2 , x 3 , X 4 , 54 (x 1 , £2 , Z3 , X 4 , «2 . P ) » Pi = U 2 
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where u 2 is the desired moment x 4 . It may be considered as the new control variable. 

(3) F 2 is the force function subject to moment trim. 

F 2 (xi,x 2 ,x i ,x 4 ,u 2 ,p) = f 2 [xi,x 2 ,x 3 ,x 4 ,g 4 (xi,x 2 ,x 3 ,x 4 ,u 2 ,p),p] ( 3 - 11 ) 


This is F* in equation (3-5). The functions g 4 and F 2 may be used to simplify the system in 
equation (3-8) as follows. 

x i = h(x\,x 2 ,p) 

x 2 F 2 (*^T > x 2 , x 4 ^ \l 2 , pPj 

X Z = X A (3-12) 

x 4 = u 2 

y = x 1 

In addition, we assume that F 2 is invertible with respect to X 3 . This partial inverse will be denoted by 

*3 = 92(xi,x 2 ,ui,x 4 ,u 2 ,p) 

F 2 [xi,x 2 ,g 2 (xi,x 2 ,ui,x 4 ,u 2 ,p),x 4 ,u 2 ,p] = m ^ 3 ' 13 ^ 

where u\ is interpreted as the desired acceleration x 2 . 


(4) The output y is x\, which will be interpreted as the generalized position of the aircraft, and x 2 
and x 2 as generalized velocity and acceleration, respectively. 


The block diagram of the model is shown in figure 2. The presence of feed-forward is an indication 
of the presence of zero dynamics. If F 2 is independent of u 2 and x 4 , then there is no feed-forward and 
no zero dynamics. This form will be referred to as pure feedback. If the model is in the pure feedback 
form, then the guidance trajectory can be obtained by inversion as follows. Suppose that the output 
y and its four time derivatives are given as functions of time, as are the parameter and its three time 
derivatives: 


Then, let 


and obtain x 


( 3 ] 

2 


by computing 


d 4l = (pPM 2) J 3> J 4) ) 

(3-14) 

pi 3 ] = (p(0),p(l), p (2) )p (3)) 

[ 4 ] [ 4 ] 

— Vc (3-15) 

x 2 = 9i(xi,x^\p) (3-16) 



Figure 2. Local model of the controlled process. 
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II Hill: 


and three of its time derivatives; then obtain by computing, with nulled-out zero dynamics, 


x 3 = g 2 (xi,X 2 , %2 


(3-17) 


and its two time derivatives, so that 


XA = x\ 

( 2 ) 

U 2 = 

U = 0 4 (xi, X 2 , X3, X 4 ,U 2 ,p) 


(3-18) 


The resulting x = (xi,x 2 ,x 3 ,x 4 ) state and control u are then taken as the guidance state and control. 


x° = x 

IT —U 


(3-19) 


The algorithm outlined above will be denoted as Tq. In practice Tq is rather complicated. The 
basic force and moment functions (//, / 6 m ) in equation (3-1) require many lines of code to implement. 
Many functions are involved— not only elementary functions such as sums, products, exponents s 
and roots, but also multivariable polynomials, vector cross products, direction cosine matrices, and 
others The functions are deeply nested in the sense of function-of-function evaluations. Furthermore, 
the transformation of (//,/ 6 m ) into (hJ 2 ,h) in equation (3-8) for a given operating mode requires 
further nonlinear computations such as the extraction of Euler angles from direction cosine matrices, 
the conversion of Cartesian coordinates into spherical, and other computations. Then it is necessary to 
obtain the inverses (^, 52 , 54 ). and, finally, pass several time derivatives through them. The coding of 
jr Q i s ver y error prone. However, the application of dynamic forms as described in reference 7 simp 1 es 
the construction and coding of T 0 from (//, f b n ) to the level at which the inverison approach becomes 
quite tractable and error free. In the remainder of this report we assume that fq is available. 

Next, suppose that zero dynamics are present. Let the acceleration error be defined as follows. 

^ = ^2(*l»*2.®3»*4.«2.P)- i 2 ( 3 ' 20) 

For the pure feedback solution the error in the command due to the ignored zero dynamics is then 

V’O = F 2 {xi,x < 2,x%x { l,v%,p) - ±2 ( 3 ‘ 21 ) 


A simple, practical way to control the effects of this error, if the zeros are weak is to close a 
loop around it by means of a regulator as shown in figure 3. The guidance regulator (Greg) provides 
stable tracking. The plant model state and control are taken as the guidance command {x ,u ), w ic 
is executable in the sense that it is a stable solution of the system (eq. (3-8)). The guidance command 
serves as input to the real plant regulator (Preg), which closes the loop on the errors between x and 
the (estimated) real plant state x?. In many practical cases this simple scheme works well. In the 
present report we are interested in cases in which the maneuvers are aggressive enough relative to the 
zero dynamics for this simple approach to be inadequate. The following example provides a simple 

illustration of such a case. 
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Figure 3. Structure of the guidance generator based on 


Linear Case 

Consider a linear case of the system in equation (3-8): 

±i = x 2 

= h = ^21^1 + ^22^2 + ^23^3 + A 24 x 4 + A 2 qu + A 2 qP 
x 3 = x 4 (3-22) 

^4 = fi = A 4 \X\ + A 42 x 2 + ^43^3 + ^44^4 + A 4 q U + A 4 qp 

V = X 1 

We assume that A 4 § is not singular, and change control coordinates from u to u 2 (see eq. (3-10)): 

u = g 4 (x, u 2 ,p) = A^ (u 2 - A 41 xi - A 42 x 2 - ^43x3 - ^442:4 - 2t 46 p) (3-23) 

Then the state equation in new coordinates becomes 


±1 = x 2 

x 2 = F 2 = C\x\ + C 2 x 2 + C3X3 + C4X4 + C5U2 + Cqp 

±3 = x 4 (3-24) 

±4 = u 2 

where 

C 5 = A 25 A^ 

c i = A 2i - C 5 A 4i (3-25) 


where 1 < 1 < 6 and i ± 5, and where u 2 , the commanded angular acceleration, is the new directly 
accessible control variable. Equation (3-24) is the linear version of equation (3-12). 

Next, assuming nonsingular C 3 , the force trim map is given by the following (u x is the commanded 

x 2)- 

x 3 = 92 = C'g’Vi - C\x\ - C 2 x 2 - C 4 x 4 - C5U2 - Cqp) (3-26) 








The pure feedback solution (x°,u§) defined by equations (3-15) through (3-19) is given by 

xj = Vc 

^2 = Vc 

x 0 = c 3 - 1( 9 < 2 ) - C 2 yP - C,yP - C 6 p(°>) (3-27) 

A = Cp(yP - C 2 y { ? - C lV P - CspO)) 

u o = c -\ y p _ c 2V P - C lV P - C 6 p< 2 >) 

The corresponding acceleration error (eq. (3-21)) is 
ipQ = C A C^ l {y { c ] - C 2 y? ) ~ C\y { c ] - C 6 p {l) ) + CsC^iyf* - C 2 yf ] - C x y { ? - C 6 p (2) ) (3-28) 

Note that ipo > s continuous for small C4 and C5, and ipo = 0 if C4 and C5 are both zero. We turn next 
to a numerical example. 


Constant Linear Example - Part I 

Suppose that the system (eq. (3-22)) is scalar, x x £ R , x G R , and 



/ 0.000 

-0.100 32.200 

0.000 

-3.220 

10.000 \ 

(x) 

U 


V 0.000 

-0.020 -1.000 

-2.000 

1.000 

-0.100; 

\p) 


In our interpretation, one radian in x 3 generates one g acceleration x 2 , and -0.1 g is generated per unit 
of control u. After the inversion (eq. (3-10)), the force function in equation (3-11) becomes 


F 2 (x,u 2 ,p) = (0.000 -0.164 28.980 


-6.440 -3.220 9.678) 


x \ 

u 2 

P 


(3-30) 


There are transmission zeros: translational acceleration x 2 is affected by angular acceleration 
command u 2 at -0.1 5/rad/sec 2 . 

We wish to transfer the system from hover at y = 0 to hover at y = 1, 000 ft in 14 sec (from t - 4 
to t = 18) as shown in figure 4, while changing the configuration parameter p from 0 to 1 and back to 
0 as shown in the figure. The complete commanded maneuver y c is composed of nine segments. Lac 
segment is a 9-degree polynomial in t with continuous splices to order 4. The pure feedback solution 

(x°,«2) and the resulting acceleration error 

■00 = F 2 (x°,u 2 ,p) - y c ( 3 ' 31 ) 
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caused by the zero dynamics are shown in figure 5. Errors of 9 ft/sec 2 are produced. We try to regulate 

the error by means of the plant model and regulator (Greg) structure shown in figure 3. A realistic 
regulator gain is 

K — { 0.0883 0.1354 6.8376 3.4005) (3-32) 

It places the two closed-loop pole pairs at 0.8 and 2 rad/sec, both with 0.5 damping. The resulting 
closed-loop error response is shown in figure 6. The acceleration error produced by the zero dynamics 
results, in this case, in a tracking error x\ that reaches 11 ft. Furthermore, the regulator amplifies the 
closed-loop acceleration error ^ • We consider next a way for reducing the tracking error by providing 
a better guidance command ( x c ,u c ). We will obtain an order-of-magnitude improvement (see fig. 13). 
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0 
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0 
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0 6 12 18 24 
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Figure 4. Commanded evolution of the output y and parameter p. 
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4 Local Zero Dynamics 


Let us return to the general case of the multiaxis system in equation (3-12). Suppose that in order 
to reduce the error ip 0 defined by equation (3-21) we introduce perturbations (771 , 772, %) along the pure 
feedback solution (a; 0 , it®). Then equation (3-20) becomes 

ip = F 2 (x*l,x%,x% + m,X4 + V2,v% + ri 3 ,p) (4-1) 

or, using equation (3-21), 

= F 2 {xi,X2,xl + 771,2:4 + 772,7^2 + *73, P) - F2(xi,X2,x®,x%U2,p) + ip 0 (4-2) 


We assume that for small perturbations, the linear part dominates, so that (approximately) 

i’ = C3771 + C4V2 + <^5% + V’O ( 4 - 3 ) 


where the Jacobian matrices 


^3 = gg 


C4 = g 


( 4 - 4 ) 




are evaluated along (x°,u%). We note that the Jacobian matrices can be computed directly from (/ 2 , / 4 ) 
in equation ( 3 - 8 ): 

C '5 = |&(^)" 1 

( 4 - 5 ) 

C i =^ r c ^ 3=3 ’ A 

If for each t in our interval of interest, T= [ti,t 2 ], we select the perturbations (77!, 772, 773) so that 

+ C^{t)r]2(t) + C§{t)rj2,{t) = - 7^0 (t) ( 4 - 6 ) 

for 0 < 7 < 1, then we have improvement at each point of T : 


iP(t) = (1 - 7)^o(f) 

The corresponding improved trajectory (rc 1 ,^) is then 


( 4 - 7 ) 




x\ = 


x 3 = ^3 + Vl 


x\ = x° A + 772 


( 4 - 8 ) 


u 2 ~ u 2 + V3 
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and the improved error (eq. (3-20)) is 


ipl = F 2 (x l ,u\,p) - ±2 


If the perturbations were free, we could have moved down in many ways, but of course the available 
perturbations are not free: according to equation (3-12) we have nonholonomic constraints, namely. 


m = m 
m = m 

so that, in fact, equation (4-6) is a differential equation, 

C 5 (t)j) + Ci(t)r) + C 3 {t)r] = - 7^(0 

where 77 = 771 . If C 5 is nonsingular, then another form of this equation is given by 

77 + A 4 (*)t) + A 3 (£)t7 = 


or, equivalently. 


*)-(-* -a 4 )G) + (b 2 )'-^ 


(4-10) 


where, of course, B 2 = C 5 l , A 4 — B 2 C±, and A 3 — B 2 C^, that is, 




Aj - 


dh d. 


Xj axj 


3 = 3,4 


We will refer to all three forms of the differential equation as the local zero dynamics (of the system 
in eq. (3-12)). We seek the stable, “steady state” solution of the local zero dynamics. We develop 
methods for the computation of such solutions in the next two sections, first using the time domain 
approach and then using the frequency domain approach. Assuming that the steady state solution can 
be computed, the proposed algorithm for the computation of the improved guidance command (x , u 2 ) 
can be summarized as follows. 

Step 1. Compute the pure feedback solution (x°, u®) and the error 7/10- 
Step 2. Compute the Jacobian matrices C 3 , C 4 , and C 5 along (x°, u®)- 
Step 3. Compute the steady state solution of the local zero dynamics. 

Step 4. Update the trajectory and compute the improved error ip. 

Step 5. Repeat steps 2 -4 until \ip\ becomes acceptable; the result is (x c ,u9>). 
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If the system (eq. (3-8)) is linear, a choice of 7 = 1 will produce the result in one iteration. 
Otherwise, 7 must be small enough for the linear part of the perturbation to be dominant. The effect of 
nonlinearity is to increase the number of required iterations. 

The differential equation (4-11) is closely related to the Isidori zero dynamics (ref. 8), which for 
the system (eq. (3-12)) are defined by the generally nonlinear differential equation 

F 2( x l ,X2>z,z,z,p) = £2 (4-15) 


b 1x1 z ^> be lts stab,e > steady state solution. The linear part of the perturbation along z(t) is given 

C5 (<)*? + C4(t)f] + Cz(t)r] = ip(t) (4-16) 

where the Jacobian matrices C* are evaluated along z(t). Thus, the homogeneous parts of equations (4-11) 
and (4-16) are the same for perturbations along z(t), but unless F 2 is linear the Jacobian matrices will 
be different in the early stages of our iteration. 

There are three approaches to solving the differential equation (4-15). 

(1) In reference 3 the approach is to represent the command y c and the parameter p as outputs 
of a relaxing exosystem, and then to design a regulator that forces the combined system to satisfy 
equation (4-15) asymptotically. 

(2) In reference 6 the approach is to solve equation (4-15) by means of a Picard-like iteration. 

(3) Our approach is to solve equation (4-15) by replacing it with a sequence of linear problems; 
each linear problem is solved by means of a Picard-like iteration (ref. 6). The advantage of this approach 
is the freedom provided by 7: it controls the intensity of the forcing function 7 ip. 

Approaches (2) and (3) are distinct only for nonlinear systems. For linear systems they are the 
same. 

The remaining problem is to develop a procedure to carry out step 3 in the above algorithm. For 
the case of our numerical example in equation (3-29), this requires finding the particular solution of the 
following differential equation. 


— 3.2207) — 6.440?) + 28.980 p — —tpQ (4-17) 

Its eigenvalues (Ai,A 2 ) = (-4.162,2.162). There is a nonminimum phase zero at 2.162. In order to 
reduce the errors shown in figures 5 and 6, we need the steady state solution of this equation. That type 
of problem is solved next. We follow the procedure given in reference 6. 
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5 Time Domain Solution of the Local Zero Dynamics 


Consider the constant linear system 


x = Ax + Bu (5-1) 

where x G R n , u G R m , and A has no eigenvalues on the ju; axis. We are not interested in the usual 
initial value problem in which x(0) and u are given and the problem is to find the solution x{t). Instead, 
we want to compute the particular solution x p for a given u. For example, what is the constant x for a 
given constant u, or what is the amplitude and phase of x for a given sinusoidal u, or what is x for an 
arbitrary given u? The answers to the first two questions are easy: 

x p(ju) = ( jujl - A)~ l Bu(joj) (5-2) 

The third question is addressed in the rest of this section and in the next section. 

Constant Linear Systems 

Assume for simplicity that A has distinct eigenvalues. Change coordinates 


to diagonalize A 


so that 


X = Pz 

(5-3) 

P~ 1 AP = A 

(5-4) 

z = Az + Tu 

(5-5) 


where T = P l B. The transition matrix for this system, which is diagonal, we separate into stable and 
unstable parts, < p~{t ) and ip + (t), respectively: 


tp(t) = e At = (p (£) + tp + (t) (5-6) 

The bilateral (ref. 6) transition matrix for equation (5-5) is then given by 

^(t) = <p~(t)U(t) - ip + (t)U(—t) (5-7) 

where U is the unit step function. The particular solution is given by the convolution 

r+OO , 

zp(t) = / <£*(£ — t)Tu(t)<It (5-8) 

J— 00 


That z p is a solution of equation (5-5) can be checked by differentiation. Note that <p (0) 4 -<p + (0) = 
I, and that and ip + cannot have nonzero entries in the same row. 

By transforming back to the natural coordinates, we obtain the bilateral transition matrix 

</>*(£) = P V ? ± (£)P“ 1 (5-9) 


p m. 
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and the bilateral input-output impulse response 

h ± (t) = P(f^{t)P~ l B (5-10) 

The particular solution to equation (5-1) is given by the convolution 

/ -fOC , 

h ± (t - t)u{t)<1t (5-11) 

-oo 

or, equivalently, the computationally more convenient form 

/ -foe i 

h ± (— a)u(t + a)do (5-12) 

-oo 


In the present report we are interested in the second-order system given by equation (4-13). 
Application of equation (5-12) with u = -7^ produces the following particular solution: 

tt l M) = -'<C h±{ -° m+<T)ia (5 - 13) 


For the scalar case, 


fj + + a ZV = l>2 • (-7"0) 


with Ai < 0 < A2 


so that 


and 



h ± (t) = 


b 2 

A2 - Ai 




( 5 - 14 ) 

( 5 - 15 ) 

( 5 - 16 ) 

(5-17) 

(5-18) 


The first component of /i ± (f), 

hf (() = + e^ l U(-t)\ 

A2 - Ai 


( 5 - 19 ) 


is shown in figure 7 for the case of equation (4-17), where (Ai,A2) = (—4.162,2.162). Note the 
rapid roll-off on either side of t = 0. One may expect that, for reasonable forcing functions —7 ip, 
the convolution may be restricted to, say, only four times the slowest time constant, which in the case 
shown is approximately 2 sec. 
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Time (sec) 

Figure 7. Numerical example of (£). 

If we define 

M l(0 = /-ooe Al(t - r) u(r)dr 

M2 (t) = /+°° e A2(t-r) M ( r ) dr 

then the particular solution of equation (5-14) is given by 

V (t) = +fi2(t)] 

V(t) = ^-[Ai/iiW+Aa^W] 
ij(t) = -a$n(t) - a 4 7)(f) - &27^(f) 


(5-20) 


(5-21) 


In order to explore this solution, suppose that we have a table of integrals so that for a given ip, 

j e Xt ip(t)dt = tf(A, f), #(A, t) = e xt ip(t ) (5-22) 


Then 


Ml(0 = e Al ^(-Ai,<) 

/42 (t) = -e A 2^(_A 2 ,f) 


(5-23) 


If ip(t) = 1, then ^(A,£) = e A /A, and 



M = J - 2 

so that the particular solution is, as expected from equation (5-14), 




~ b 27 


(5-24) 


77 = 0 (5-25) 

m = o 
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If ip(t) = sin(aif), then 


so that 
where 


H\ = j[-Ai sin(wt) - ucos(ut)] 

u 2 = -e A2 *tf(-A 2 ,£) = 5 [-% sin (ut) - wcos(wi)] 

' ^2 tW 


77 : — 7asin(u;f + 6 ) 


ae = 


^2 


-w 2 + joja .4 + 03 


(5-26) 


(5-27) 

(5-28) 


Discontinuities in ip or its derivatives spawn leading (noncausal) and trailing transients. For exam- 


ple, consider a pulse, 


( 0, t <t\ 

lf)(t) = < lpo(t), t\ < t < f 2 

l o, t 2 <t 


(5-29) 


Then 


( 0, t < t\ 

ji\ = J e Al *[^ 0 (-Ai,f) - $o(-Ai,ti)], h<t<t 2 
[ e Al *[^o(-Ai,t 2 ) - ^o(-Ai,fi)], f 2 < t 


(5-30) 


and 




' e A2< ['Fo(-A 2 ,f 2 ) - ®o(-A 2 ,*i)]> 1 < h 
< e A2i [^o(-A 2 ,f 2 ) - A 2 , i)], £i < t < t 2 

,0, f 2 < t 


(5-31) 


That is, there is a pair of transients at either end of the pulse. One member of the pair leads the 
discontinuity, the other trails it. For example, suppose that 


( 0, t < 0 

ip{i) — J 1 — cos(u>£), 0 < t < t 2 = (5-32) 

l 0, f 2 < t 

One way to get the particular solution is to splice three segments of particular solutions. 


r 0, t < 0 

rj(t ) = < oq — ocos(u;£ + 0), 0 < t < t 2 (5-33) 

1 0, t 2 <t 

where a 0 = -f>27/ a 3> anti a and ^ are §' ven b y equation (5-28). But that simple approach produces 
undesirable discontinuities in 77 and r) at t = 0 and t — t 2 ■ On the other hand, the particular solution 


t l(t) = + f* 2 (t)} (5-34) 

A 2 - Ai 

and its derivative are everywhere continuous. The leading and trailing transients satisfy the homogeneous 
part 

i ] + 0477 + 0377 = 0 (5-35) 

These transients prepare the system state x and control u for the approaching discontinuity, and the 
transients are invisible at the system output y. Thus, the bilateral solution has definite advantages. 
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Numerical Solution 


Numerical solutions of the convolution integrals are needed in practical applications. The integrals 
in equation (5-20) may be computed approximately as follows. 

Ml = e Ali jlooe-^rjdr = E£=o I^nT-T e~ X ' T ^dr « 

Ml = — E£L 0 (e AlT H^ - nT - T) + rftt - nT )]/ 2 

(5-36) 

M2 = e A2 * Jt 00 e~ X2T ip{r)dT = E£L 0 ti+nT+T e~ X2T 'tp{r)dr « 

M 2 = j^^^TAT) EjLo(e~ A2T ) n [^ + nT + T) + ^ + nT)]/2 

The approximations have been normalized to give correct results for ip = 1. Consequently, the approx- 
imate solution to our problem is given by 

il = + /4) 

^2 = Jj-A, OlMl + ^2^2) 

^3 = -OSHf - fl 4»72 “ 6 27^(^) 

For example, consider yet again the system discussed previously, 

(03,04,6) = (-9.000 2.000 0.311) 

(Ai,A 2 ) = (-4.162 2.162) 


(5-37) 


(5-38) 


Assume that the forcing function is 


(—7 ip) = 20sin(u;£ -I- 7t/4) 


so that the exact solution is easily obtained as the particular solution: 


vU u ) = 


-0.311 

(— + Juj2 — 9) 




(5-39) 


(5-40) 


We take the sampling period T = 0.05 sec and provide a look-ahead/look-back of N = 40 samples 
(±2 sec). The approximate computation 77* (0) is compared to the exact 77(0) for several frequencies. A 
frequency scan from [0, 10] rad/sec is shown figure 8. The approximate and exact solutions are practically 
indistinguishable in the figure. Hence this type of algorithm is effective for constant systems. We now 
turn to the time-varying case. 
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Figure 8. Exact and approximate solutions. 


Linear Time- Varying Systems 

In general, the local zero dynamics given by equation (4-13) may have time-varying coefficients. 
Therefore, consider the following model 

x = [Aq + 6A(t)\x + B(t)u(t) (5-41) 

where x e R n , u € R m , and 8 A and B are some reasonable functions of time. The iteration in refer- 
ence 6, restricted to linear systems, will converge to the particular solution if the following conditions 
hold: 


(1) ^4q has no eigenvalues on the ju axis; 0 ± (f) is the bilateral solution for the system 


x = Aqx 


and let the norm of be 

. roo . 

11^*11 = J_ 00 \^fk^\ dt 


(2) There is a constant K so that for all — oo < t < oo 


max |6ajo(f)| < K 
i,j 


(5-42) 

(5-43) 


(5-44) 


(3) The product 

kuH < i 


(5-45) 
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The solution is constructed by means of a Picard-like iteration: 


x 0 (t) = 0 

®l(0 = f™oo - r)S(r)u(r)dr 

x m+ 1(0 = xi(t) + /~o ^(t - r)tfi4(r)x m (r)dr 


The difference at each iteration 


so that 


e m (f) = x m+ i(<) - x m (f) = f^oo 0 ± (f - r)M(r)e m _ 1 (r)(ir 
max t ||e m (t)|| < max t ||xi(O||(i ; f||0 ± ||) m 


x p (t) 


lim 

m— >oo 


x m{t) 


(5-47) 


(5-48) 


In actual applications, it may be advantageous to change coordinates. For example, suppose that 
Aq has only real and distinct eigenvalues, and that P diagonalizes it. Then 

11^11= (5-49) 

that is, the sum of time constants, and condition (3) becomes 

<1 (5-50) 


In the scalar case of equation (4- 1 2), separate the time variation as follows 

V + [«04 + <5a 4 (f)]77 + [<303 + Sa 3 (t)]r] = -b 2 (t)'y'ip(t) (5-51) 

If the eigenvalues of Aq satisfy Ai < 0 < A 2 , and if for —00 < t < 00 

max[|6a 3 (f)|, |<5a 4 (f)|](ri + r 2 ) < 1 (5-52) 

then the iteration will converge to the particular solution. 

In summary, if the coefficients of the local zero dynamics vary in time sufficiently little relative to 
their average value, then the iteration will produce the particular solution, and so step 3 in our algorithm 
can be carried out. The overall algorithm is thus complete. The time domain approach seems to require 
a lot of bookkeeping with respect to the eigenvalue pattern. The frequency domain approach discussed 
next is simpler, especially for multiaxes systems. 


25 




6 Frequency Domain Solution of the Local Zero Dynamics 

Consider again the local zero dynamics (eq. (4-11)) 

C 5 (t)rj + C 4 (t)f] + C 3 (t)r] = -7 ip(t) (6-1) 

We are interested in obtaining the particular solution. In the preceding section, we took the time domain 
approach, which has some drawbacks. The need to distinguish between root patterns is a nuisance. 
Also, we would like to consider cases where C5 drifts through zero, which leads to changing order 
of the differential equation. The frequency domain approach, discussed next, is in some ways much 
simpler to implement. 

Let us rewrite equation (6-1) as a sum of constant and time-varying linear operators 

(L 0 + 6L)r] = -ytj) (6-2) 

where D = d/dt and 

Lq = C$D 2 + C$D + Cj 

(6-3) 

6L = 6C 5 {t)D 2 + 6C 4 (t)D + 6C s {t) 

Let 7 ] 0 be the solution of the time invariant part 

Vo = Lq 1 (6-4) 

and change coordinates so that 77 = 770 + ?7i . Then 

(L 0 + 8L)(t]q + rn) = -y*l) (6-5) 

That is, 

(L 0 + 6L)r]i = -<5Lt / 0 (6-6) 

This is the old equation (6-2) with a new forcing term, which is the error caused by tjq . After i 
iterations, the correction is 

r) j+ i = -Lq'ILti, (6-7) 

If ||io : 'Mil < 1, then the iterations 

V = J2 Vi ( 6 - 8 ) 

converge to the solution of equation (6-2) and hence of the original equation (6-1). Consider next the 
construction of Lq 1 . 

The Fourier transform of equation (6-4) is 

VO O’w) = G 0 {juj)[-'rip(j(j)] (6-9) 

where the transfer function 

Go{ju>) = [C®(jlo) 2 + C 4 ju + C3] 1 (6-10) 
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The forcing function i>{t) is available only at sampling points, t}){ nT ). The Shannon reconstruction 
(ref. 9), with sinc(x) = sin(x)/x. 


+ OO 

1>(t) = T £ nT)sinc[uj s {t - nT)/2]/T (6-11) 

— OO 

converts the sequence ip(nT) into an analog signal, ip(t), that is band-limited to |w| ^ u>g/2 tt/T. 
The Fourier transform of sinc{u s t/2)/T is 



1, M < w s f 2 
0, jcuj > oJs/2 


Hence, the transform of the output of G 0 with input sinc(u s t/2)/T is 


= G 0 {juj)S(ju) 


( 6 - 12 ) 


(6-13) 


and its inverse is 

or, equivalently, because of symmetries, 


(6-14) 


%(t) = 


1 
7 r 



cos ut — Goy sin ut)dw 


(6-15) 


where Gq x and Gq v are, respectively, the real and imaginary parts of Gq. The integral is then computed 
numerically (we used an IMSL, Inc., routine (ref. 10)) for 2 N s + 1 values of t to obtain the sequence 


{g 0 (kT), —N s < k < N s } 


(6-16) 


Then, finally, the frequency domain approximation ff 0 to the particular solution of equation (6 4) is 
given by the following discrete convolution: 


N s 

r f{nT) = -1T £ iP[(n-k)T)}g 0 (kT) (6-17) 

k=-N s 


The derivatives if 1 and t) u are computed similarly by means of equation (6-17), except that go 
is replaced by g\ and g<i, which are the inverse transforms of G\(juj) = juGo(ju) and = 

(ju>) 2 Go(juj), respectively. 

The potential aliasing problem arising from the multiplication in forming 8Lru may be avoided by 
means of a noncausal filter Gj , which provides roll-off without phase lag. In that case the combined 

transfer function 

Go(ju)) = [C'fO'cj) 2 + C®ju> + C3] Gf(ju) (6-18) 

where Gf(s) is symmetric about both real and imaginary axes. For example, if, as before, (c 5 , c 4 , c 3 ) = 
(-3.22, -6.44, 28.98), we may choose 

Go(jto) = [-3.22 (ju) 2 - 6A4(ju) + 28.98]- 1 [l + WS) 4 ]- 1 (6-19) 
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The response g Q of the (unstable) G 0 to sinc{u> s t/2)/T is shown in figure 9. Also shown is 
the corresponding bilateral impulse response hf (see fig. 7). It may be noted that the two curves are 
practically indistinguishable, which is not surprising since hf(t) and G 0 (juj) (without the filter) are a 
Fourier transform pair. 



Time (sec) 

Figure 9. Comparison of hf with the response go of Gq to sinc(w s t/2)/T. 


Thus we have two effective ways for computing the steady state solution of the local zero dynamics. 
One solves the problem in time domain producing the approximation rf\ the other solves the problem 
in frequency domain producing the approximation rf° . 
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7 Outline of the Proposed Algorithm 


An approach has been presented in this report for the computation of the guidance commands 
(x c , u c ) that guide the system along y c . The system may be multiaxis, time- varying, and nonlinear, but 
of the following type (see eq. (3-8)) 


±1 fl (^T ) p) 

±2 = /2(®li®2»*3.®4.«»P) 

±3=0:4 (7-1) 

±4 = h{xi,X2,Xz,X4,U,p) 

y = x 1 


where x^,u E R m and p E Rr is a time-dependent parameter. This is a realistic model of an aircraft. 
The algorithm consists of the following steps. 

Step 1. Change the control variable to u 2 = fi{x\,X 2 ,x^,X 4 ,u,p) so that the model becomes 


±1 = /l(®l,*2*P) 

±2 = F 2 (xi,X2,X3,X4,u 2 ,p) 

±3 = £4 (7-2) 

±4 = u 2 

y = xi 


Step 2. Initialize (z = 0) the following iteration by inverting the pure feedback part of the 

system. 


x\ = y c 

= 9l(x\,x\,p) 

= g 2 (x\,xl,x l 2 ,0,0,p) 


x 4 ~ x 3 
— x\ 


Step 3. Compute the error in x 2 caused by the neglected dynamics: 


(7-3) 


= F 2 {x\ , x l 2 , x\, x\ , u\ , p) ~ ±2 (7-4) 

Skip Steps 4-6 if |^| is small enough. 

Step 4. Compute the Jacobian matrices C3, C\, and C5 along [x i ^u^). 

Step 5. Compute iteratively the particular solution of the local zero dynamics, 

C\fi + C\fj + C\p = —7 ipi (7-5) 
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Step 6. Update trajectory (x\u l 2 ) 


x 

x 

x 

X 

u 


i+1 _0 

J - x l 
i+1 _ T 0 
2 ~ X 2 

3 +1 = x 3 + V 

\ +1 =X\ + T) 

2 +1 = u 2 + j i 


Step 7. Compute the control in natural coordinates 


u l = 94 (A ,x l 2 ,x l 3 ,x\, u \ , p) 


(7-6) 


(7-7) 


The result is the guidance trajectory ( x c , u c ), which includes the effects of zero dynamics. For 
linear systems one iteration with 7 = 1 will suffice. For nonlinear systems, 7 is chosen small enough to 
ensure the dominance of the linear approximation. If many iterations are needed in a given problem, then 
the accumulation of small errors may make the iterate (x\,x\,u l 2 ) inconsistent with the nonholonomic 
constraints, 

r ui _ Ji 

H - x 4 (7-8) 

x\ = u\ 

That can be easily rectified by passing each iterate through a noncausal observer for the constraint 
(eq- (7-8)). 

The resulting structure of the command generator shown in figure 10 is the same as in figure 3, 
except that the pure feedback guidance (x°,u°), which does not account for zero dynamics, is replaced 
by (x c , u c ), which does. We now apply, by way of illustration, the complete algorithm to compute the 
guidance command ( x c ,u c ) for several examples. 


pPl.yJ 4 ! 



Disturb 


yP 


Figure 10. Structure of the guidance generator based on ( x c ,u c ). 
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8 Examples 


In this section we apply the numerical solution of the local zero condition to the linear model 
discussed previously and describe the resulting closed-loop behavior. Later in this section, we apply the 
proposed guidance method to a time-varying system, and then to a nonlinear system. 

Constant Linear Example - Part II 

In the first part of this example (see Constant Linear Example — Part I) we have computed and 
stored in memory the error t[) o(£) on the time interval [0,26] at 20 samples per sec. Now we pass this 
sequence ipo(nT) with 7 = 1 through the algorithm (eqs. (5-36) and (5-37)), 81 samples (N = 40) at a 
time at each sample, to obtain the corrections ( 77 *,?)*, 7 /*). The first panel of figure 11 shows both time 
domain 77 * and frequency domain rf (with N s = 40) solutions. They are indistinguishable. The next 
two panels show the approximations to 77 and 77 , where 

fi**(n) = [ri*(n + l)-v*(n-l))/(2T) 

7 f» = [ 77 * (n + 1 ) - 77 *(n - 1 )]/( 2 T) (8_1) 


The pairs of curves in the bottom panels are again practically indistinguishable; hence we have 
consistency: ( 77 *, 77 *) are good approximations to the time derivatives of ( 77 *, 77 *). It may be noted that 
there are leading and trailing transients, such as before t = 4 and after t = 18. That effect is even 
more pronounced in figure 12 , where improved command x c (solid line) is compared with the original 
command x° that ignores the zero dynamics. In the figure the pure feedback solution is shown dotted. 
It is a copy from figures 5 and 6 . A noteworthy beneficial effect is that the corrected command is 
smaller, smoother, and gentler than x°. The tracking errors are shown in figure 13. They have been 
reduced by a factor of 20. Thus, we obtain better tracking by means of gentler commands. 

.3 
0 
-.3 


n* n *J\ 


- A 




J 





Figure 1 1 . Numeric solution of the zero dynamics. 


33 










Linear Time- Varying Example 


The state equation in this example is 


0 

1 

0 

0 ^ 


( 0 \ 

0 

-0.1 

32.2 

-10.18 

X + 

Pit) 

0 

0 

0 

1 

0 

\0 

0 

0 

0 ) 


l 1 / 


The corresponding zero dynamics are given by (see eq. (4-11)) 

p(t)i) — 10.18?) + 32.277 = —ipo 


( 8 - 2 ) 


(8-3) 


The root locus for —3.22 < p < 3.22 is shown in figure 14. At p — —3.22 there is a saddle, one 
root being close to —5 and the other close to +2. As p increases, the stable root goes off to —00 (at 
p = 0 there is a drop in order!) and reappears from +00 while the unstable root passes through 3. The 
roots meet for p = 0.8, resulting in a repeated unstable pair near +6. Then the zero dynamics turn into 
an unstable spiral. At p = 3.22, the roots rest at (1.58 ±j'2.74). 
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-10 

P 

= -3.22 

p = 0 

j 







-10 -5 0 5 10 

Figure 14. Root locus. 


We wish to guide the system along the same trajectory y c as shown in figure 4, except that the 
parameter p is changed as shown in figure 15. The system error response without correcting for zero 
dynamics is shown in figure 16, and the system error response with corrections for zero dynamics is 
shown in figure 17. It should be noted that an order of magnitude improvement is obtained despite time 
variation of the zero dynamics. The corrections were implemented in frequency domain. The following 
noncausal transfer function was used to compute ??(f) w . 
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Figure 15. Variation of p. 
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G 0 (juj) = [p(juf - 10.18(jw) + 32.2]“ 1 [1 + (iu/8) 4 ] -1 (8-4) 

where p is locally frozen (see below). The derivative was computed by means of G\(ju>) — 

jujGo(juj). Rapid frequency roll-off without phase lag was provided by the noncausal filter with the 
pole pattern shown in figure 18. The filter requires a look-ahead/look-back of 0.5 sec, but since we are 
already prepared to look-ahead/look-back 2 sec, the noncausal nature of the filter is not a problem. The 
g 0 (kT),-AO < k < 40, needed in the convolution for r/ 61 ' 1 , and gi(kT), —40 < k < 40, needed for rf ^ , 
were computed every fifth sample using the value of p at the beginning of each quintet. 



Figure 1 8. Roll-off filter poles. 


Nonlinear Time- Varying Example 

The model used next is motivated by the fact that generally the aircraft force and moment generators 
stiffen with speed. We multiply the linear force and moment used in the preceding example by a factor 
that is quadratic in speed x 2 to obtain the following nonlinear state equation. 


±1 = X2 

x 2 = h = i a 2l x l + a 22 x 2 + «23^3 + <*24^4 + «25« + «26P)[1 + («27 3; 2) 2 ] 

±3 = X4 (8-5) 

±4 = h = (<*41^1 + a 42^2 + a 43 x 3 + a 44 X 4 + a 4 5 U + a467>)[l + («47^2) 2 ] 

y = x 1 


The maneuver to be executed is the same as before (fig. 4). With <127 = 047 = 0.01, the multiplying 
factor ranges from 1 to 2.5. Figure 19 shows the spectrum (Ai, A2) of the zero dynamics with Jacobian 
matrices frozen at each sample. Both zeros move out with increasing X 2 ■ Also shown are the regulator 
gains, varied so that the closed-loop poles of the perturbation model with frozen Jacobian matrices at 
each sample remain, as before, at 0.8 and 2 rad/sec, both with 0.5 damping. The tracking errors after 
only one iteration (7 = 1) are shown in figure 20; they are small. It is particularly noteworthy that 
the acceleration error of the corrected guidance ip (solid line) is an order of magnitude smaller than the 
pure feedback error ipQ (dotted line). 
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9 Conclusion 


An algorithm for computing state trajectories and controls that guide a nonlinear system through 
given sequences of control points has been presented. It is assumed that the control point schedule is 
known several points in advance. In the algorithm, the control points are linked by polynomial segments. 
The first approximation of the guiding state trajectory is obtained by equating the system output with 
the polynomial schedule and inverting the pure feedback part of the system. In the presence of zero 
dynamics this first approximation will produce acceleration errors. A correction is obtained by means 
of the stable, steady state solution of the local zero dynamics. If nonminimum phase zeros are present, 
future values of the approximate solution are used. Numerical tests indicate that an order-of-magnitude 
reduction in the acceleration error is possible with look-ahead of only four times the time constant of 
the unstable zero. Therefore, the control points must be available that much in advance. Usually, this 
does not pose a problem in the fully automatic mode of operation of the system, since the higher-level 
planner that provides the control point schedule has a much wider horizon than the servosystem being 
discussed here. Although a regulation process must always be causal, a guidance process may employ 
noncausal computations. 
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